// 'ob' == offset binary
// 'tag' == Write line numbers before each line output
//
//  TODO: 目前的版本中，从bin文件提取“波形”数据、CDS处理、写入fits文件等功能，全部在
//          同一个函数中实现，代码显得臃肿，不便于后期的更新维护。


#include "ccd290_extract.h"

void extract_ccd_data_ob_tag( char fname[], char ofits[], ccd290_config* cfg ) {

    printf("--> reading binary data from: %s\n", fname);

    unsigned short udata_E, udata_E_;
    unsigned short udata_F, udata_F_;
    unsigned short udata_G, udata_G_;
    unsigned short udata_H, udata_H_;

//    int16_t data_E, data_F, data_G, data_H;

    FILE *fb = fopen(fname,"rb");

    int num = cfg->LEN*4;	// total number of chars
    unsigned short *regs = (unsigned short*)malloc(sizeof(unsigned short)*num);
    int num_chk = fread(regs, sizeof(unsigned short), num, fb);

    if( num_chk != num ){
    	printf("Error in reading binary data: %s\n", fname);
    }

    fclose(fb);

    int32_t *ch_E, *ch_F, *ch_G, *ch_H;
    ch_E = (int32_t*)malloc(sizeof(int32_t)*cfg->LEN);
    ch_F = (int32_t*)malloc(sizeof(int32_t)*cfg->LEN);
    ch_G = (int32_t*)malloc(sizeof(int32_t)*cfg->LEN);
    ch_H = (int32_t*)malloc(sizeof(int32_t)*cfg->LEN);

    int i,j;
    int cnt=0;

#ifdef _PRINT_A_LINE_
    int imin = (_DEBUG_LINE_NUM_*(cfg->N_COL + cfg->HAS_TAG) * cfg->SAMP_NUM_PER_REG*4);    // for debug
#endif

#ifdef _WRITE_LINE_TAG_
    FILE *fp_tag = fopen("line_tag.txt","w");
#endif
    
	int is_tag;
    for(i=0; i<num; i+=4) {

        udata_E = regs[i+0];
        udata_F = regs[i+1];
        udata_G = regs[i+2];
        udata_H = regs[i+3];
        
		is_tag = ( i%((cfg->N_COL + cfg->HAS_TAG) * cfg->SAMP_NUM_PER_REG * 4) ) < cfg->SAMP_NUM_PER_REG*4;

        j = i+4;
        if( (j < num) ) {

        //  ==================
        //  重要！！！
		//	需要注意EFGH的顺序
        //  ==================
	        uint16_t E, F, G, H;
			if( is_tag == 1 ) {
		        ch_E[cnt] = udata_E;
		        ch_F[cnt] = udata_F;
		        ch_G[cnt] = udata_H;
		        ch_H[cnt] = udata_G;

#ifdef _WRITE_LINE_TAG_
		        fprintf(fp_tag,"%10d %10d %10d %10d\n", udata_E, udata_F, udata_G, udata_H);
#endif
			} else {
		        udata_E_ = regs[j+0];
		        udata_F_ = regs[j+1];
		        udata_G_ = regs[j+2];
		        udata_H_ = regs[j+3];

		    //	shift 4-bits to recover original LTC2271 AD outputs
	        //  unsigned short E, F, G, H;
		        E = ((udata_E & 0x0fff) << 4) | ((udata_E_ & 0xf000) >> 12);
		        F = ((udata_F & 0x0fff) << 4) | ((udata_F_ & 0xf000) >> 12);
		        G = ((udata_G & 0x0fff) << 4) | ((udata_G_ & 0xf000) >> 12);
		        H = ((udata_H & 0x0fff) << 4) | ((udata_H_ & 0xf000) >> 12);

		        // data_E = fun2s(E);
		        // data_F = fun2s(F);
		        // data_G = fun2s(G);
		        // data_H = fun2s(H);

		        // ch_E[cnt] = data_E;
		        // ch_F[cnt] = data_F;
		        // ch_G[cnt] = data_H;
		        // ch_H[cnt] = data_G;

                ch_E[cnt] = E;		//实际对应G通道
		        ch_F[cnt] = F;		//实际对应H通道
		        
		    //	注意通道顺序：G和H需要交换
		        ch_G[cnt] = H;		//实际对应F通道
		        ch_H[cnt] = G;		//实际对应E通道

			}

#ifdef _PRINT_A_LINE_
			E = udata_E;
			F = udata_F;
			G = udata_G;
			H = udata_H;
            if( i>=imin && i<imin+(cfg->N_COL + cfg->HAS_TAG) * cfg->SAMP_NUM_PER_REG*4){
            	printf(" %6d %6d %6d %6d 0x%04x 0x%04x 0x%04x 0x%04x\n", data_E, data_F, data_G, data_H, E, F, G, H);
                // printf("*%6u %6u %6u %6u 0x%04x 0x%04x 0x%04x 0x%04x\n", data_E, data_F, data_G, data_H, E, F, G, H);
            }
#endif
            cnt++;
        }
    }

#ifdef _WRITE_LINE_TAG_
    fclose(fp_tag);
#endif
    
//	CDS
    uint16_t *ccd = (uint16_t*)malloc(sizeof(uint16_t)*cfg->N_ROW*(cfg->N_COL + cfg->HAS_TAG)*4);
    
    for(i=0; i<cfg->N_ROW*(cfg->N_COL+cfg->HAS_TAG)*4; i++) 
        ccd[i] = 0;

//  get settings from ccd290_config * cfg
    // double size1 = cfg->size1;
    // double size2 = cfg->size2;
	uint idx_offset = cfg->idx_offset; //-5;

    printf("idx_offset = %d\n", cfg->idx_offset);

    uint E_idx1 = cfg->E_idx1;
    uint E_idx2 = cfg->E_idx2;
    uint E_idx3 = cfg->E_idx3;
    uint E_idx4 = cfg->E_idx4;

    uint F_idx1 = cfg->F_idx1;
    uint F_idx2 = cfg->F_idx2;
    uint F_idx3 = cfg->F_idx3;
    uint F_idx4 = cfg->F_idx4;
    
    uint G_idx1 = cfg->G_idx1;
    uint G_idx2 = cfg->G_idx2;
    uint G_idx3 = cfg->G_idx3;
    uint G_idx4 = cfg->G_idx4;
    
    uint H_idx1 = cfg->H_idx1;
    uint H_idx2 = cfg->H_idx2;
    uint H_idx3 = cfg->H_idx3;
    uint H_idx4 = cfg->H_idx4;
    
    printf("E_idx1 = %2d, E_idx2 = %2d, E_idx3 = %2d, E_idx4 = %2d\n", E_idx1, E_idx2, E_idx3, E_idx4);
    printf("F_idx1 = %2d, F_idx2 = %2d, F_idx3 = %2d, F_idx4 = %2d\n", F_idx1, F_idx2, F_idx3, F_idx4);
    printf("G_idx1 = %2d, G_idx2 = %2d, G_idx3 = %2d, G_idx4 = %2d\n", G_idx1, G_idx2, G_idx3, G_idx4);
    printf("H_idx1 = %2d, H_idx2 = %2d, H_idx3 = %2d, H_idx4 = %2d\n", H_idx1, H_idx2, H_idx3, H_idx4);
  
    int E_ref = 0, E_sig = 0, E_both = 0;
    int F_ref = 0, F_sig = 0, F_both = 0;
    int G_ref = 0, G_sig = 0, G_both = 0;
    int H_ref = 0, H_sig = 0, H_both = 0;

    printf("--> running CDS processing\n");
    
    int row_cnt=0, col_cnt=0;
    int idx;

	/************************************************************************
		CCD数据矩阵对应的通道分布，为了使图像拼接完整，在下面的代码中需要交换G、H通道的数据
		E | F
		--|--
		H | G
	*************************************************************************/

    for(idx=idx_offset; idx<cfg->LEN; idx+=cfg->SAMP_NUM_PER_REG) {
    
	    int idxE = row_cnt*(cfg->N_COL + cfg->HAS_TAG) * 2 + col_cnt;
	    int idxF = row_cnt*(cfg->N_COL + cfg->HAS_TAG) * 2 + ((cfg->N_COL + cfg->HAS_TAG) * 2 - col_cnt - 1);
	    int idxG = (cfg->N_ROW*2 - row_cnt - 1)*(cfg->N_COL + cfg->HAS_TAG) * 2 + col_cnt;
	    int idxH = (cfg->N_ROW*2 - row_cnt - 1)*(cfg->N_COL + cfg->HAS_TAG) * 2 + ((cfg->N_COL + cfg->HAS_TAG) * 2 - col_cnt - 1);

		    
	    int sumE1=0, sumE2=0;
	    int sumF1=0, sumF2=0;
	    int sumG1=0, sumG2=0;
	    int sumH1=0, sumH2=0;

        int E1_ok = 1, E2_ok = 1;
        int F1_ok = 1, F2_ok = 1;
        int G1_ok = 1, G2_ok = 1;
        int H1_ok = 1, H2_ok = 1;

		int cnt;

		int PIX_MAX = 65535;
		int PIX_MIN = 1;

		for( cnt=E_idx1; cnt<=E_idx2; cnt++ ){
            int tmp = (int)ch_E[idx+cnt];
            if( tmp >= PIX_MAX )
                E1_ok = 0;
            
            sumE1 += tmp;
        }
            
		for( cnt=E_idx3; cnt<=E_idx4; cnt++ ){
            int tmp = (int)ch_E[idx+cnt];
            if( tmp <= PIX_MIN )
                E2_ok = 0;

            sumE2 += tmp;
        }

		for( cnt=F_idx1; cnt<=F_idx2; cnt++ ){
            int tmp = (int)ch_F[idx+cnt];
            if( tmp >= PIX_MAX )
                F1_ok = 0;
            sumF1 += tmp;
        }

		for( cnt=F_idx3; cnt<=F_idx4; cnt++ ) {
            int tmp = (int)ch_F[idx+cnt];
            if( tmp <= PIX_MIN )
                F2_ok = 0;
            sumF2 += tmp;
        }

		for( cnt=G_idx1; cnt<=G_idx2; cnt++ ) {
            int tmp = (int)ch_G[idx+cnt];
            if( tmp >= PIX_MAX )
                G1_ok = 0;
            sumG1 += tmp;
        }

		for( cnt=G_idx3; cnt<=G_idx4; cnt++ ) {
            int tmp = (int)ch_G[idx+cnt];
            if( tmp <= PIX_MIN )
                G2_ok = 0;
            sumG2 += tmp;
        }

		for( cnt=H_idx1; cnt<=H_idx2; cnt++ ) {
            int tmp = (int)ch_H[idx+cnt];
            if( tmp >= PIX_MAX )
                H1_ok = 0;
            sumH1 += tmp;
        }

		for( cnt=H_idx3; cnt<=H_idx4; cnt++ ) {
            int tmp = (int)ch_H[idx+cnt];
            if( tmp <= PIX_MIN )
                H2_ok = 0;
            sumH2 += tmp;
        }

		/*TODO:
			add offset to ccd[idx#] if ( sum#1/size1 - sum#2/size2 ) is negative
		*/
	    int32_t ccd_E = ( ( sumE1/cfg->E_size_ref - sumE2/cfg->E_size_sig ) + cfg->cds_offset );
	    int32_t ccd_F = ( ( sumF1/cfg->F_size_ref - sumF2/cfg->F_size_sig ) + cfg->cds_offset );
	    int32_t ccd_G = ( ( sumG1/cfg->G_size_ref - sumG2/cfg->G_size_sig ) + cfg->cds_offset );
	    int32_t ccd_H = ( ( sumH1/cfg->H_size_ref - sumH2/cfg->H_size_sig ) + cfg->cds_offset );
		
    //  assign values according to X1_ok & X2_ok flags

		if( (E1_ok==1) & (E2_ok==1) ){
			ccd[idxE] = (uint16_t)ccd_E;
		} else if( (E1_ok==0) && (E2_ok==1) ) {
			ccd[idxE] = 65534;
            E_ref++;
		} else if( (E1_ok==1) && (E2_ok==0) ) {
            ccd[idxE] = 32767;
            E_sig++;
        } else if( (E1_ok==0) && (E2_ok==0) ){
            ccd[idxE] = 65535;
            E_both++;
        }

		if( (F1_ok==1) & (F2_ok==1) ){
			ccd[idxF] = (uint16_t)ccd_F;
		} else if( (F1_ok==0) && (F2_ok==1) ) {
			ccd[idxF] = 65534;
            F_ref++;
		} else if( (F1_ok==1) && (F2_ok==0) ) {
            ccd[idxF] = 32767;
            F_sig++;
        } else if( (F1_ok==0) && (F2_ok==0) ){
            ccd[idxF] = 65535;
            F_both++;
        }

		if( (G1_ok==1) & (G2_ok==1) ){
			ccd[idxG] = (uint16_t)ccd_G;
		} else if( (G1_ok==0) && (G2_ok==1) ) {
			ccd[idxG] = 65534;
            G_ref++;
		} else if( (G1_ok==1) && (G2_ok==0) ) {
            ccd[idxG] = 32767;
            G_sig++;
        } else if( (G1_ok==0) && (G2_ok==0) ){
            ccd[idxG] = 65535;
            G_both++;
        }

		if( (H1_ok==1) & (H2_ok==1) ){
			ccd[idxH] = (uint16_t)ccd_H;
		} else if( (H1_ok==0) && (H2_ok==1) ) {
			ccd[idxH] = 65534;
            H_ref++;
		} else if( (H1_ok==1) && (H2_ok==0) ) {
            ccd[idxH] = 32767;
            H_sig++;
        } else if( (H1_ok==0) && (H2_ok==0) ){
            ccd[idxH] = 65535;
            H_both++;
        }
		
    	if( (cfg->HAS_TAG == 1) && (col_cnt == 0) ){
		    ccd[idxE] = ch_E[idx];
		    ccd[idxF] = ch_F[idx];
		    ccd[idxG] = ch_G[idx];
		    ccd[idxH] = ch_H[idx];
    	}
        
    	col_cnt++;  // 必须先+1,然后再进行判断
        if( col_cnt >= (cfg->N_COL + cfg->HAS_TAG) ) {
            col_cnt -= (cfg->N_COL + cfg->HAS_TAG);
            row_cnt++;
        }
    }

    printf("******************** CDS error statistics ********************\n");
    printf("* X_ref:  ref level exceeds 32767, sig level is OK\n");
    printf("* X_sig:  sig level is lower than -32767\n");
    printf("* X_both: ref level exceeds 32767 & sig level is lower than -32767\n");

    printf("=> E_ref = %8d,", E_ref);
    printf(" E_sig = %8d,", E_sig);
    printf(" E_both = %8d\n", E_both);

    printf("=> F_ref = %8d,", F_ref);
    printf(" F_sig = %8d,", F_sig);
    printf(" F_both = %8d\n", F_both);

    printf("=> G_ref = %8d,", G_ref);
    printf(" G_sig = %8d,", G_sig);
    printf(" G_both = %8d\n", G_both);

    printf("=> H_ref = %8d,", H_ref);
    printf(" H_sig = %8d,", H_sig);
    printf(" H_both = %8d\n", H_both);

//	===================================================
//  write CDS process data into fits
    fitsfile *fptr;
    int status;
    long fpixel = 1, naxis = 2, nelements, exposure;
    
    long naxes[2] = {(cfg->N_COL + cfg->HAS_TAG)*2, cfg->N_ROW * 2};
    
    if( cfg->include_line_num == 0 ){
    	naxes[0] = cfg->N_COL * 2;
    	naxes[1] = cfg->N_ROW * 2;
    }

    status = 0;
    
	fits_create_file(&fptr, ofits, &status); /* create new file */

    /*Create the primary array image (16-bit short uinteger pixels) */
    /*	fits_create_img(fptr, SHORT_IMG, naxis, naxes, &status);*/
    fits_create_img(fptr, USHORT_IMG, naxis, naxes, &status);

    /*Write a keyword; must pass the ADDRESS of the value */
    exposure = 3.;
    fits_update_key(fptr, TLONG, "EXPOSURE", &exposure, "Total Exposure Time", &status);

    long speed = 636;
    fits_update_key(fptr, TLONG, "RD_SPD", &speed, "READOUT SPEED (KHz)", &status);

    double temp = -96.5;
    fits_update_key(fptr, TDOUBLE, "TEMP", &temp, "CCD operating temperature", &status);

    double v_rst_gate = 12;
    fits_update_key(fptr, TDOUBLE, "RST_GATE", &v_rst_gate, "RESET GATE (V)", &status);

    double v_rst_dr = 16.25;
    fits_update_key(fptr, TDOUBLE, "RST_DR", &v_rst_dr, "RESET DRAIN (V)", &status);

    double v_out_drain = 30.5;
    fits_update_key(fptr, TDOUBLE, "OUT_DRN", &v_out_drain, "OUTPUT DRAIN (V)", &status);
    
    double v_dp_drain = 29;
    fits_update_key(fptr, TDOUBLE, "DP_DRAIN", &v_dp_drain, "DUMP DRAIN (V)", &status);

    double v_op_1v = 0.5;
    fits_update_key(fptr, TDOUBLE, "OP_1V", &v_op_1v, "OP_1V (V)", &status);


    // 写入头文件信息
    


    /* Write the array of integers to the image */
    nelements = naxes[0] * naxes[1];
    
    if( cfg->include_line_num == 0 ) {
		uint16_t *ccd_tmp = (uint16_t*)malloc(sizeof(uint16_t) * cfg->N_ROW * cfg->N_COL * 4);
		int m,n;
		for(m=0; m<2*cfg->N_ROW; m++){
			for(n=1; n<=2*cfg->N_COL; n++){
				// ccd_tmp[m*2*cfg->N_COL+n] = ccd[m*2*(cfg->N_COL+cfg->HAS_TAG) + n];
                ccd_tmp[m*2*cfg->N_COL+ (n-1)] = ccd[m*2*(cfg->N_COL+cfg->HAS_TAG) + n]; // ccd_tmp[]中的n忘记-1了
			}
		}
		
		fits_write_img(fptr, TUSHORT, fpixel, nelements, ccd_tmp, &status);
		free(ccd_tmp);
	} else {
	    fits_write_img(fptr, TUSHORT, fpixel, nelements, ccd, &status);
	}

    fits_close_file(fptr, &status);

    printf("--> data saved into fits: %s\n", ofits);

    fits_report_error(stderr, status);
//	===================================================

    free(ch_E);
    free(ch_F);
    free(ch_G);
    free(ch_H);

    free(regs);
    free(ccd);
}


